      SUBROUTINE DIAGL(U,V,IPAR)
      IMPLICIT NONE
*     .. Parameters ..
      INTEGER LDA
      PARAMETER (LDA=1)
      INTEGER LOCLEN
      PARAMETER (LOCLEN=LDA)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*),V(*)
      INTEGER IPAR(*)
*     ..
*     .. Arrays in Common ..
      DOUBLE PRECISION Q1(LOCLEN)
*     ..
*     .. External Subroutines ..
      EXTERNAL DCOPY,DVPROD
*     ..
*     .. Common blocks ..
      COMMON /PIMQ1/Q1
*     ..
      CALL DCOPY(IPAR(4),U,1,V,1)
      CALL DVPROD(IPAR(4),Q1,1,V,1)
      RETURN

      END
      SUBROUTINE DIAGR(U,V,IPAR)
      IMPLICIT NONE
*     .. Parameters ..
      INTEGER LDA
      PARAMETER (LDA=1)
      INTEGER LOCLEN
      PARAMETER (LOCLEN=LDA)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*),V(*)
      INTEGER IPAR(*)
*     ..
*     .. Arrays in Common ..
      DOUBLE PRECISION Q2(LOCLEN)
*     ..
*     .. External Subroutines ..
      EXTERNAL DCOPY,DVPROD
*     ..
*     .. Common blocks ..
      COMMON /PIMQ2/Q2
*     ..
      CALL DCOPY(IPAR(4),U,1,V,1)
      CALL DVPROD(IPAR(4),Q2,1,V,1)
      RETURN

      END
      SUBROUTINE PRINTV(N,U)
      IMPLICIT NONE
*     .. Scalar Arguments ..
      INTEGER N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*)
*     ..
*     .. Local Scalars ..
      INTEGER I
*     ..
      WRITE (6,FMT=9000) (U(I),I=1,N)
      RETURN

 9000 FORMAT (4(D14.8,1X))
      END
      SUBROUTINE PROGRESS(LOCLEN,ITNO,NORMRES,X,RES,TRUERES)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      DOUBLE PRECISION NORMRES
      INTEGER ITNO,LOCLEN
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION RES(*),TRUERES(*),X(*)
*     ..
*     .. External Subroutines ..
      EXTERNAL PRINTV
*     ..
*     WRITE (6,FMT=9000) ITNO,NORMRES
*     WRITE (6,FMT=9010) 'X:'
*     CALL PRINTV(LOCLEN,X)
*     WRITE (6,FMT=9010) 'RES:'
*     CALL PRINTV(LOCLEN,RES)
*     WRITE (6,FMT=9010) 'TRUE RES:'
*     CALL PRINTV(LOCLEN,TRUERES)
      RETURN
 9000 FORMAT (/,I5,1X,D16.10)
 9010 FORMAT (/,A)
      END
      SUBROUTINE REPORT(NAME,IPAR,DPAR,ET,X)
      IMPLICIT NONE
*     .. Scalar Arguments ..
      REAL ET
      CHARACTER*(*) NAME
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION DPAR(*),X(*)
      INTEGER IPAR(*)
*     ..
*     .. External Subroutines ..
      EXTERNAL PRINTV
*     ..
      WRITE (6,FMT=9000) NAME,IPAR(2),IPAR(11),IPAR(12),IPAR(13),
     +  DPAR(2),ET
      CALL PRINTV(IPAR(4),X)


      RETURN

 9000 FORMAT (/,A,': N=',I6,/,2X,'k=',I6,' status=',I6,' steperr=',I6,
     +       ' norm=',D16.10,/,2X,'Execution time(s)= ',F15.8)
      END
      SUBROUTINE DULSCOEF(M,C)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      INTEGER M
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION C(*)
*     ..
      IF (M.EQ.1) THEN
          C(1) = 7.0D0/6.0D0
          C(2) = 5.0D0/6.0D0

      ELSE IF (M.EQ.2) THEN
          C(1) = 35.0D0/32.0D0
          C(2) = 25.0D0/16.0D0
          C(3) = 35.0D0/32.0D0

      ELSE IF (M.EQ.3) THEN
          C(1) = 37.0D0/40.0D0
          C(2) = 49.0D0/40.0D0
          C(3) = 91.0D0/40.0D0
          C(4) = 63.0D0/40.0D0

      ELSE IF (M.EQ.4) THEN
          C(1) = 91.0D0/96.0D0
          C(2) = 7.0D0/12.0D0
          C(3) = 21.0D0/16.0D0
          C(4) = 7.0D0/2.0D0
          C(5) = 77.0D0/32.0D0

      ELSE IF (M.EQ.5) THEN
          C(1) = 117.0D0/112.0D0
          C(2) = 87.0D0/112.0D0
          C(3) = -3.0D0/7.0D0
          C(4) = 33.0D0/28.0D0
          C(5) = 627.0D0/112.0D0
          C(6) = 429.0D0/112.0D0

      ELSE IF (M.EQ.6) THEN
          C(1) = 1059.0D0/1024.0D0
          C(2) = 687.0D0/512.0D0
          C(3) = 429.0D0/1024.0D0
          C(4) = -759.0D0/256.0D0
          C(5) = 429.0D0/1024.0D0
          C(6) = 4719.0D0/512.0D0
          C(7) = 6435.0D0/1024.0D0

      ELSE IF (M.EQ.7) THEN
          C(1) = 1117.0D0/1152.0D0
          C(2) = 1397.0D0/1152.0D0
          C(3) = 979.0D0/384.0D0
          C(4) = -143.0D0/1152.0D0
          C(5) = -10153.0D0/1152.0D0
          C(6) = -715.0D0/384.0D0
          C(7) = 17875.0D0/1152.0D0
          C(8) = 12155.0D0/1152.0D0

      ELSE IF (M.EQ.8) THEN
          C(1) = 2497.0D0/2560.0D0
          C(2) = 451.0D0/640.0D0
          C(3) = 143.0D0/80.0D0
          C(4) = 4147.0D0/640.0D0
          C(5) = -143.0D0/256.0D0
          C(6) = -13871.0D0/640.0D0
          C(7) = -2431.0D0/320.0D0
          C(8) = 17017.0D0/640.0D0
          C(9) = 46189.0D0/2560.0D0

      ELSE IF (M.EQ.9) THEN
          C(1) = 2879.0D0/2816.0D0
          C(2) = 2249.0D0/2816.0D0
          C(3) = -923.0D0/1408.0D0
          C(4) = 4537.0D0/1408.0D0
          C(5) = 6253.0D0/352.0D0
          C(6) = 221.0D0/704.0D0
          C(7) = -69173.0D0/1408.0D0
          C(8) = -29393.0D0/1408.0D0
          C(9) = 130169.0D0/2816.0D0
          C(10) = 88179.0D0/2816.0D0

      ELSE IF (M.EQ.10) THEN
          C(1) = 4173.0D0/4096.0D0
          C(2) = 2587.0D0/2048.0D0
          C(3) = 169.0D0/4096.0D0
          C(4) = -3107.0D0/512.0D0
          C(5) = 12597.0D0/2048.0D0
          C(6) = 48841.0D0/1024.0D0
          C(7) = 12597.0D0/2048.0D0
          C(8) = -54587.0D0/512.0D0
          C(9) = -205751.0D0/4096.0D0
          C(10) = 499681.0D0/6144.0D0
          C(11) = 676039.0D0/12288.0D0

      ELSE IF (M.EQ.11) THEN
          C(1) = 13081.0D0/13312.0D0
          C(2) = 15853.0D0/13312.0D0
          C(3) = 36643.0D0/13312.0D0
          C(4) = -32657.0D0/13312.0D0
          C(5) = -163591.0D0/6656.0D0
          C(6) = 72029.0D0/6656.0D0
          C(7) = 818159.0D0/6656.0D0
          C(8) = 178619.0D0/6656.0D0
          C(9) = -3000347.0D0/13312.0D0
          C(10) = -1508087.0D0/13312.0D0
          C(11) = 1924111.0D0/13312.0D0
          C(12) = 1300075.0D0/13312.0D0

      ELSE IF (M.EQ.12) THEN
          C(1) = 28243.0D0/28672.0D0
          C(2) = 1363.0D0/1792.0D0
          C(3) = 30209.0D0/14336.0D0
          C(4) = 34901.0D0/3584.0D0
          C(5) = -267767.0D0/28672.0D0
          C(6) = -73321.0D0/896.0D0
          C(7) = 15181.0D0/1024.0D0
          C(8) = 546193.0D0/1792.0D0
          C(9) = 2503573.0D0/28672.0D0
          C(10) = -839477.0D0/1792.0D0
          C(11) = -3528775.0D0/14336.0D0
          C(12) = 928625.0D0/3584.0D0
          C(13) = 5014575.0D0/28672.0D0
      END IF

      RETURN

      END
      SUBROUTINE DVPROD(N,DX,INCX,DY,INCY)
      IMPLICIT NONE
*
*     Modified from daxpy level 1 BLAS
*     element-wise vector multiplication, y<-x*y
*     Rudnei Dias da Cunha, 16/6/93
*
*     constant times a vector plus a vector.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*
*     .. Scalar Arguments ..
      INTEGER INCX,INCY,N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION DX(*),DY(*)
*     ..
*     .. Local Scalars ..
      INTEGER I,IX,IY,M,MP1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC MOD
*     ..
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
          DY(IY) = DY(IY)*DX(IX)
          IX = IX + INCX
          IY = IY + INCY
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD(N,4)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          DY(I) = DY(I)*DX(I)
   30 CONTINUE
      IF (N.LT.4) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
          DY(I) = DY(I)*DX(I)
          DY(I+1) = DY(I+1)*DX(I+1)
          DY(I+2) = DY(I+2)*DX(I+2)
          DY(I+3) = DY(I+3)*DX(I+3)
   50 CONTINUE
      RETURN

      END
      SUBROUTINE DWLSCOEF(M,C)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      INTEGER M
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION C(*)
*     ..
      IF (M.EQ.1) THEN
          C(1) = 6.0D0/5.0D0
          C(2) = 4.0D0/5.0D0

      ELSE IF (M.EQ.2) THEN
          C(1) = 8.0D0/7.0D0
          C(2) = 12.0D0/7.0D0
          C(3) = 8.0D0/7.0D0

      ELSE IF (M.EQ.3) THEN
          C(1) = 8.0D0/9.0D0
          C(2) = 4.0D0/3.0D0
          C(3) = 8.0D0/3.0D0
          C(4) = 16.0D0/9.0D0

      ELSE IF (M.EQ.4) THEN
          C(1) = 10.0D0/11.0D0
          C(2) = 4.0D0/11.0D0
          C(3) = 16.0D0/11.0D0
          C(4) = 48.0D0/11.0D0
          C(5) = 32.0D0/11.0D0

      ELSE IF (M.EQ.5) THEN
          C(1) = 14.0D0/13.0D0
          C(2) = 8.0D0/13.0D0
          C(3) = -16.0D0/13.0D0
          C(4) = 16.0D0/13.0D0
          C(5) = 96.0D0/13.0D0
          C(6) = 64.0D0/13.0D0

      ELSE IF (M.EQ.6) THEN
          C(1) = 16.0D0/15.0D0
          C(2) = 8.0D0/5.0D0
          C(3) = 0.0D0
          C(4) = -16.0D0/3.0D0
          C(6) = 0.0D0
          C(7) = 64.0D0/5.0D0
          C(8) = 128.0D0/15.0D0

      ELSE IF (M.EQ.7) THEN
          C(1) = 16.0D0/17.0D0
          C(2) = 24.0D0/17.0D0
          C(3) = 64.0D0/17.0D0
          C(4) = -16.0D0/17.0D0
          C(5) = -256.0D0/17.0D0
          C(6) = -64.0D0/17.0D0
          C(7) = 384.0D0/17.0D0
          C(8) = 256.0D0/17.0D0

      ELSE IF (M.EQ.8) THEN
          C(1) = 18.0D0/19.0D0
          C(2) = 8.0D0/19.0D0
          C(3) = 48.0D0/19.0D0
          C(4) = 208.0D0/19.0D0
          C(5) = -32.0D0/19.0D0
          C(6) = -704.0D0/19.0D0
          C(7) = -256.0D0/19.0D0
          C(8) = 768.0D0/19.0D0
          C(9) = 512.0D0/19.0D0

      ELSE IF (M.EQ.9) THEN
          C(1) = 22.0D0/21.0D0
          C(2) = 4.0D0/7.0D0
          C(3) = -16.0D0/7.0D0
          C(4) = 16.0D0/3.0D0
          C(5) = 32.0D0
          C(6) = 0.0D0
          C(7) = -256.0D0/3.0D0
          C(8) = -256.0D0/7.0D0
          C(9) = 512.0D0/7.0D0
          C(10) = 1024.0D0/21.0D0

      ELSE IF (M.EQ.10) THEN
          C(1) = 24.0D0/23.0D0
          C(2) = 36.0D0/23.0D0
          C(3) = -24.0D0/23.0D0
          C(4) = -304.0D0/23.0D0
          C(5) = 256.0D0/23.0D0
          C(6) = 2048.0D0/23.0D0
          C(7) = 256.0D0/23.0D0
          C(8) = -4352.0D0/23.0D0
          C(9) = -2048.0D0/23.0D0
          C(10) = 3072.0D0/23.0D0
          C(11) = 2048.0D0/23.0D0

      ELSE IF (M.EQ.11) THEN
          C(1) = 24.0D0/25.0D0
          C(2) = 36.0D0/25.0D0
          C(3) = 24.0D0/5.0D0
          C(4) = -32.0D0/5.0D0
          C(5) = -256.0D0/5.0D0
          C(6) = 512.0D0/25.0D0
          C(7) = 5888.0D0/25.0D0
          C(8) = 256.0D0/5.0D0
          C(9) = -2048.0D0/5.0D0
          C(10) = -1024.0D0/5.0D0
          C(11) = 6144.0D0/25.0D0
          C(12) = 4096.0D0/25.0D0

      ELSE IF (M.EQ.12) THEN
          C(1) = 26.0D0/27.0D0
          C(2) = 4.0D0/9.0D0
          C(3) = 32.0D0/9.0D0
          C(4) = 544.0D0/27.0D0
          C(5) = -64.0D0/3.0D0
          C(6) = -512.0D0/3.0D0
          C(7) = 256.0D0/9.0D0
          C(8) = 1792.0D0/3.0D0
          C(9) = 512.0D0/3.0D0
          C(10) = -23552.0D0/27.0D0
          C(11) = -4096.0D0/9.0D0
          C(12) = 4096.0D0/9.0D0
          C(13) = 8192.0D0/27.0D0
      END IF

      RETURN

      END
      SUBROUTINE MATVEC(U,V,IPAR)
      IMPLICIT NONE
*     .. Parameters ..
      INTEGER LDA
      PARAMETER (LDA=1)
      INTEGER LOCLEN
      PARAMETER (LOCLEN=LDA)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*),V(*)
      INTEGER IPAR(*)
*     ..
*     .. Arrays in Common ..
      DOUBLE PRECISION MatrizSist(LDA,LOCLEN)
*     ..
*     .. External Subroutines ..
      EXTERNAL DGEMV
*     ..
*     .. Common blocks ..
      COMMON /PIMA/MatrizSist
*     ..
      CALL DGEMV('N',IPAR(2),IPAR(2),1.0D0,MatrizSist,IPAR(1),U,
     1           1,0.0D0,V,1)
      RETURN

      END


      SUBROUTINE PDSUM(ISIZE,X,IPAR)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
*     .. Scalar Arguments ..
      INTEGER ISIZE
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
      RETURN

      END


      DOUBLE PRECISION FUNCTION PDNRM2(LOCLEN,U,IPAR)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
*     .. Scalar Arguments ..
      INTEGER LOCLEN
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL DNRM2
*     ..
      PDNRM2 = DNRM2(LOCLEN,U,1)
      RETURN

      END


      SUBROUTINE POLYL(U,V,IPAR)
      IMPLICIT NONE
*     .. Parameters ..
      INTEGER LDA
      PARAMETER (LDA=1)
      INTEGER LOCLEN
      PARAMETER (LOCLEN=LDA)
      INTEGER NGAMMA
      PARAMETER (NGAMMA=13)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*),V(*)
      INTEGER IPAR(*)
*     ..
*     .. Scalars in Common ..
      INTEGER MMM
*     ..
*     .. Arrays in Common ..
      DOUBLE PRECISION GAMMA(NGAMMA),Q1(LOCLEN)
*     ..
*     .. Local Scalars ..
      INTEGER I,LOCALN
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION T(LOCLEN),W(LOCLEN)
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DSCAL,DVPROD,MATVEC
*     ..
*     .. Common blocks ..
      COMMON /B0001/GAMMA,MMM
      COMMON /PIMQ1/Q1
*     ..
      LOCALN = IPAR(4)
      CALL DCOPY(LOCALN,U,1,T,1)
      CALL DVPROD(LOCALN,Q1,1,T,1)
      CALL DCOPY(LOCALN,T,1,V,1)
      CALL DSCAL(LOCALN,GAMMA(MMM+1),V,1)
      DO 10 I = 1,MMM
          CALL MATVEC(V,W,IPAR)
          CALL DVPROD(LOCALN,Q1,1,W,1)
          CALL DAXPY(LOCALN,-1.0D0,W,1,V,1)
          CALL DAXPY(LOCALN,GAMMA(MMM-I+1),T,1,V,1)
   10 CONTINUE
      RETURN

      END

	SUBROUTINE POLYR(U,V,IPAR)
      IMPLICIT NONE
*     .. Parameters ..
      INTEGER LDA
      PARAMETER (LDA=1)
      INTEGER LOCLEN
      PARAMETER (LOCLEN=LDA)
      INTEGER NGAMMA
      PARAMETER (NGAMMA=13)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*),V(*)
      INTEGER IPAR(*)
*     ..
*     .. Scalars in Common ..
      INTEGER MMM
*     ..
*     .. Arrays in Common ..
      DOUBLE PRECISION GAMMA(NGAMMA),Q2(LOCLEN)
*     ..
*     .. Local Scalars ..
      INTEGER I,LOCALN
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION T(LOCLEN),W(LOCLEN)
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DSCAL,DVPROD,MATVEC
*     ..
*     .. Common blocks ..
      COMMON /B0001/GAMMA,MMM
      COMMON /PIMQ2/Q2
*     ..
      LOCALN = IPAR(4)
      CALL DCOPY(LOCALN,U,1,T,1)
      CALL DVPROD(LOCALN,Q2,1,T,1)
      CALL DCOPY(LOCALN,T,1,V,1)
      CALL DSCAL(LOCALN,GAMMA(MMM+1),V,1)
      DO 10 I = 1,MMM
          CALL MATVEC(V,W,IPAR)
          CALL DVPROD(LOCALN,Q2,1,W,1)
          CALL DAXPY(LOCALN,-1.0D0,W,1,V,1)
          CALL DAXPY(LOCALN,GAMMA(MMM-I+1),T,1,V,1)
   10 CONTINUE
      RETURN

      END
      SUBROUTINE TMATVEC(U,V,IPAR)
      IMPLICIT NONE
*     .. Parameters ..
      INTEGER LDA
      PARAMETER (LDA=1)
      INTEGER LOCLEN
      PARAMETER (LOCLEN=LDA)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION U(*),V(*)
      INTEGER IPAR(*)
*     ..
*     .. Arrays in Common ..
      DOUBLE PRECISION MatrizSist(LDA,LOCLEN)
*     ..
*     .. External Subroutines ..
      EXTERNAL DGEMV
*     ..
*     .. Common blocks ..
      COMMON /PIMA/MatrizSist
*     ..
      CALL DGEMV('T',IPAR(2),IPAR(2),1.0D0,MatrizSist,IPAR(1),U,
     1           1,0.0D0,V,1)
      RETURN

      END
